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Abstract 

A method to evaluate the nonrelativistic electron-repulsion, nuclear attraction and kinetic energy 
three-electron integrals over Slater orbitals appearing in Hylleraas-CI (Hy-CI) electron structure 
calculations on atoms is shown. It consists on the direct integration over the interelectronic coordi- 
nate 7y and the sucessive integration over the coordinates of one of the electrons. All the integrals 
are expressed as linear combinations of basic two-electron integrals. These last are solved in terms 
of auxiliary two-electron integrals which are easy to compute and have high accuracy. The use 
of auxiliary three-electron ones is avoided, with great saving of storage memory. Therefore this 
method can be used for Hy-CI calculations on atoms with number of electrons iV > 5. It has been 
possible to calculate the kinetic energy also in terms of basic two-electron integrals by using the 
Hamiltonian in Hylleraas coordinates, for this purpose some mathematical aspects like derivatives 
of the spherical harmonics with respect to the polar angles and recursion relations are treated and 
some new relations are given. 
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1. Introduction 

The Hylleraas- Configuration Interaction (Hy-CI) wave functions [21 13] are of great interest 
in Quantum Chemistry because they lead to highly accurate energy values, or benchmark 
calculations, which usually serve to test or calibrate other Quantum Mechanical methods. 
This possibility of achieving high accuracy has been demonstrated with recent calculations 
of the helium atom [3], the H2 molecule [5], the lithium [6], and beryllium atoms 011]. The 
application of the Hy-CI method can be extended to the calculation of the states of the second 
period of elements. Hy-CI calculations on the B atom [U] are in progress. Also Hy-CI have 
been successfully used for the calculation excited states. One example of their applications 
is the study of atomic excitations during /3~-decay processes of atoms and ions, and the 
calculation of the final state probabilities [10]. In the future Hy-CI wave functions may 
be applied to modern problems of Physical Chemistry as confined atoms [IT] , and Nuclear 
Physics like the study of excitations during nuclear reactions [12]. Some of these nuclear 
reactions are important in Nuclear Medicine, like the Boron Neutron Capture Therapy [T3] . 
a promising therapy against head and many types of cancer. For all these reasons it is 
important to provide to future workers in the field with different analytical methods to solve 
all occurring Hy-CI integrals and to compare their effectiveness in practical calculations. 

The kinetic energy integrals, which are generated by the kinetic energy operator part 
of the Hamiltonian, are not as complex as the repulsion four-electron integrals, but their 
evaluation is also difficult. Note that the two-electron kinetic-energy integrals are needed 
not only in computer programs for two-electron systems, but in any program code for larger 
atoms. 

In paper I [13] of this series we have derived analytical expressions for the three-electron 
kinetic-energy integrals occurring in the Hy-CI method. The Hamiltonian in Hylleraas 
coordinates was used |15J . The method of derivation consisted in the direct application of 
the differentiating operators on the wave function. This included to perform derivatives over 
spherical harmonics and the use of recursion relations over these functions. We distinguished 
between the cases of one quantum number m = 0, m > 0, and m < 0. In this paper we 
obtain general expressions for any m quantum number. The most important fact is that the 
method employed here can be extended in the same fashion to the three-electron case. 

The advantage of our method is due to the direct integration over the interelectronic 
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distances and the coordinates of one of the electrons, reducing the integral to a lower order. 
The three- and two-electron kinetic energy integrals are reduced to a linear combination of 
usual two-electron integrals. In appendix A the two-electron integrals are derived in detail. 

In the method by Sims and Hagstrom [19] the interelectronic distance is expanded into 
one-electron distances. The expansion of the interelectronic distances is a concept they use 
systematically in the evaluation of all kind of integrals. This fact leads to the appearance of 
three-electron auxiliary integrals W. For contrary using our method for the three-electron 
kinetic energy integrals no three-electron auxiliary integrals are needed, with a great saving 
of computer memory. Note that the auxiliary W integrals are three-fold and have to be 
calculated for a large number of powers and exponents. The number of exponents grows 
with the atomic number. 

We have also checked the two-electron kinetic energy integrals using the Kolos and 
Roothaan transformation [16], obtaining completely agreement of the integrals to more than 
30 decimal digits. Along this work we have used quadruple precision (30 decimal digits in 
our computer). In appendix B, the Kolos and Roothaan transformation for two-electron 
kinetic energy integrals is derived. If a similar transformation for the three-electron case 
can be obtained is still open. The Kolos and Roothaan transformation has shown to be a 
factor of times computationally faster than the integrals shown here, with identical memory 
requirements. 

The resulting formulas have been programed and used in highly accurate calculations 
on the ground and excited states of He atom with an accuracy > 1 • 10~ 10 a.u., to obtain 
some benchmark calculations for the , n = 2 — 4 of the ion Li + [lOj and in preliminary 
calculations on the B atom [9]. We shall see that the Kolos and Roothaan transformation 
provides an effective computational method for the two-electron case. Whereas the integrals 
and method showed here shall help us to better understanding of the underlying relations 
in the kinetic energy integrals and it can be used as a method for the three-electron case. 

2. The Hylleraas-CI wave function 

The Hy-CI wave function [21 [3] is a linear combination of configurations built up with 
Slater-type atomic orbitals with higher angular momentum as in the regular Configuration- 
Interaction (CI) procedure and configurations including one interelectronic distance as for 
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Hylleraas-type trial wave functions [T7] . The Hy-CI and CI wave functions for an n-electron 
systems are defined as: 

N 

* = £)C^, § q = d{L 2 )A<P qX , (1) 

9=1 

where $ g are symmetry adapted configurations, N is the number of configurations and the 
constants C q are determined variationally. The operator 0(L 2 ) projects over the proper 
spatial space, so that every configuration is eigenfunction of the square of the angular mo- 
mentum operator L 2 . A is the n-particle antisymmetrization operator, and \ is the spin 
eigenfunction: 

x = [(a/3 - /3a)...(a/3 - /3a)a] (2) 

where for even electron systems the last a spin function is omitted. The spatial part of the 
basis functions are Hartree products of Slater orbitals: 

n 
k=l 

The powers v take the values 0, or 1. For v = the wave function reduces effectively to 
a CI wave function. To understand the meaning of the Hy-CI wave function, note that 
even powers of the interelectronic coordinate are equivalent to products of p-, d-, ■ ■ ■ type 
one-electron orbitals, for instance: 

rlj=p{i)p{j)- (4) 

Furthermore, higher odd powers can be expressed as r^rf™. It can be demonstrated using 
the addition theorem of spherical harmonics that a configuration times is equivalent to 
an infinite expansion of higher angular momentum one-electron orbitals: 

s(i)s0>« = s Ws(j) + P(i)p{j) + d{i)d{j) + /(i)/0') + • • • ( 5 ) 

In this way a single Hy-CI configuration is equivalent to infinite sum of CI configurations 
containing the excitations of the involved electrons to all higher momentum orbitals. 

In this section the two-electron kinetic energy integrals occurring in the Hy-CI method 
when using the Hamiltonian in Hylleraas coordinates [15] shall be derived. For any atomic 
number N > 2 two kinds of kinetic energy integrals occur: two- and three-electron ones. 
The three-electron kinetic energy integrals can be found in Refs. [HI IE]- In this work 
we will perform the derivatives directly over and using Hamiltonian written in polar 
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and interelectronic coordinates, we shall prove its correctness and will solve the resulting 
integrals in terms of radial two-electron integrals. 

Let us define the Slater orbitals specified by the quantum numbers n, m and I with an 
unnormalized radial part and orthonormal spherical harmonics: 



1 w 



n'-l -a'rym' 



(6) 



The spherical harmonics in Condon and Shortley phases [20, p. 52] are given by: 

"2/ + l(Z-m)P 1/2 



Y, 



i ^> ' 



P, m (cos0)e im * 



(7) 



4tt (l + m)\ 

with the associated Legendre functions P z m (cos#). The spherical harmonics and associated 
Legendre functions used along this work are written explicitly in [211 P- 14], and defined as 
in Ref. [? ]. They obey the condition: 



;-ir*r 



(8) 



We define the one-electron charge distributions by expanding or linearizing the products 
of spherical harmonics with equal argument using the formula [? , Eq. (12)]: 



)Y l V 



i+v 

E 

L=\l-V\ 



2L + 1 



-,1/2 



C L (l',m';l,m)Y™'- m (e, 



(9) 



where the Condon-Shortley coefficients [201 Eqs. (6-11)] are defined by: 

47T 



C L (l',m',l,m) 



•ym —m ^ 



)Y 



I K u >' 



sm 



(10) 



_2L + 1_ 

Li satisfies the triangular condition \h — < L < ij + l\ and the restriction Lj > |Mj|. 
The summation is done in steps of two: Li = \U — \k — l'^ + 2, . . . , ij + l[ — 2, k + l[ and 
Mi = m'i — rrii. The lowest value of Li depends also on m^, for simplicity we will use the 
notation \l — V\ to recall on I and I'. For discussion about the use of the Condon and Shortley 
coefficients, see Appendix B and and Appendix in Ref. [? ]. 
The charge distributions are: 

i+v 

n N , L , M (r) = = E ^L + l) 1/2 C L (l',m';l,m)r N - 1 e^ r Yi I (9,<f ) ), (11) 



L=\l-l'\ 
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where N — n + n' — 1, and the exponents u = a + a'. In the next we will use capital letters 
N, L, M for the quantum numbers of charge distributions, while low letters n, I, m for the 
quantum numbers of the orbitals will be used. 
The rotation: 

Consider the triangle formed by r± and r 2 , see Figure 1. The original idea from Calais 
and Lowdin consisted on making a rotation of one coordinate axis, which allows to make 
a change of variable in the integral. Letting pass for a moment the z-axis through the 
r\ coordinate, the variables are transformed as 9 2 — > 9 12 , and <fi 2 — > (fin ■ This may be 
understood graphically in Figure 1. The volume element of electron 2 may be then written: 

dr 2 = r\dr 2 sin 6i 2 d9 12 d(fi± 2 . (12) 

As 9 12 is related to r 12 through the cosine theorem: 

r i2 — r i + r 2 ~ 2rir 2 cos#i2, (13) 

we can differentiate on the left hand side of the equation with respect to r 12 , and in the 
other side with respect to 9\ 2 : 

2r 12 dr 12 = 2r 1 r 2 sin9i 2 d9 12 , (14) 
obtaining the relation, which will be used to change the variable of integration: 

sin 9\ 2 d9\ 2 = —^dri 2 . (15) 
nr 2 

This equation will be used in the direct integration over the variable r± 2 . With the change 
of integration variable 9 12 — > r 12 , the integration domain changes as J* — > Krl^rlr There 
are two regions of integration: when r\ < r 2 then \r\ — r 2 \ = r 2 — r\ and when r 2 < r\ then 
\f\ — t 2 \ — r\ — r 2 . The separation of domains of integration proposed by Perkins [? ] has 
proved to be the best way. Other strategies lead to the same type of formulas. The domain 
of integration is divided into two parts: 

J = D 1 -D 2 (16) 

with: 

rri pri+r-2 poo i>ri+r 2 poo rr2+r\ 

Di-> dr 2 / dr 12 + dr 2 dr 12 ->■ / dr 2 dr 12 , (17) 

Jo J r\—T2 J r\ J t\— T2 Jo J r\— V2 
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roo i>r\+T2 poo rr±+r2 

D 2 ^ dr 2 / dr 12 - dr 2 dr 12 . (18) 

J n Jr\— T2 J T\ J r 2 —r\ 

The direct integration over the r 12 coordinate [? ] leads to the separation of the variables 
ri and r 2 : 

r r 1+ r 2} 1 L("+2)/2j , . 

■Ai-m 2 12 („ + l) ^ V^-V 1 ^ ' 

where \_(v + 2)/2j is the floor function and denotes the integral part of (5 + 2)/2. 

The rotation of the axis produces a rotation of the functions, in this case, spherical 
harmonics. They suffer a transformation given by the rotation matrices, which can be 
written: 




For a n-electron system the effective Hamiltonian in Hylleraas coordinates may be written 



2 ^ dr 2 ' ri dr* ^ n ^ dr 2 - ^ ru dra ' n,- 

i=l 4 i=l i=l i<j L J i<j J J i<j J 



1 ^ r i + r ij r ] & 2 1 r ij + T ik r ]k 9 2 

2^ TiVij dridrij 2^^-f ryn & dr i:j dr ik 



1^ 1 9' 1A 1 d l 1 ^ cot 9i d 

~ 2^^M ~ 2 2^ r? sin 2 Oidtf ~ 2^ rf <9ft 

1 = 1 1 1 2=1 1 J=l ' 

Z—/ \ r- r - -- sin Hi 9 



r^r^ sin 9i 2 rf / dftdr^- 



Er 7 - sin#,- ... 9 
— 7T sin (u>< - </?,) - - . (21) 



As the Hy-CI wave function consists on only one r*y per configuration and treating now 
a two electron atom/ion, the terms of the Hamiltonian including ^2- and a a \ — vanish. 
The angular momentum operator can be extracted: 



Id 2 1 A 1 d 2 1 A cot 9i d 

i 

and its eigenvalue equation used: 



V 1 f 2 = 1 V" 1 1 V" 1 1 cot 9j d 

^ rf * 2 ^ rf <9# 2 2 ^ rf sin 2 ft <90 2 2 ^ rf «9ft ' 1 j 

i=l 4 i=l 1 1 i=l 1 1 i=l 1 
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L 2 fa = h(h + 

with li the angular quantum number of the orbital fa. 

From the variational principle one obtains the matrix eigenvalue problem: 



(23) 



(H - EA)C = 0, (24) 

where the matrix elements are: 



H, 



ki 



J QtHfydr, (25) 



® k $idT. (26) 



The Hamiltonian in Hylleraas coordinates can be separated in a sum of one-electron 
potential and kinetic energy operators. 

n 

H = J2 H{1) = f (1) + V(l) = f R (l) + f e (l) + V(l). (27) 

i=i 

The kinetic energy operator can be again separated into radial and angular parts. Let us 
evaluate the kinetic energy of electron 1 and the interelectronic distance r^: 

f(l)=f fll (l)+f i?2 (l)+f i?3 (l)+f L (l)+f ei (l)+f e2 (l)+f (? 3(l). (28) 
The radial operators are: 



Id 2 Id 

= "2*T n (29) 

f Rl (l) = J ±t^zd / , (30) 

^a(l) = (31) 
r 12 or 12 

and the angular ones: 



8 



2 r{ 89 f 89i sm 9 1 8(f) j 

T 6l 1 = 2 — —JL , 33 

rxr 12 sin t^i 89 1 8r 12 

= -^cotg / 12 2 ri r2 / , (34) 
2 rfri2 89idri2 

Te 3 (±) = — ^ sin (0i - 2 ) — — — , (35) 

rinzsmdi 8(j)\8ri2 

the potential energy operator is: 

V(l) = -^ (36) 

In the two-electron kinetic energy integrals we may encounter CI configurations v = and 
configurations of Hy-CI type with v — 1. The three-electron integrals J are defined [T4"] : 



T/ AT TIT TIT -, •, \ ,?7l', 777 2 ,771, ,77l 3 ,771^ 

J{Nx, N%, N3; ui, U2, ui3] 1, l) hA / 2W ? 

are defined in Appendix A, Eq. (A. 3) and their value can be calculated using Eqs. 
(A.12,A.27). 

The nuclear attraction potential energy is: 



I PE {1) = (^(r!)^^)^)^^!^^)^^)^^)^} 

ry 7/ jit tit at -i 1 N 7771 ,777',, 7772 .TTlo, 7713 , 7710 /Qrr\ 

= -ZJ(A^ 1 ,^ 2 ,Ar 3 ;w 1 ,w 2 ,a;3;l,l) il/iii2 1 4>/3ii / ! (37) 

with Z the atomic nuclear charge. 
The radial kinetic energy integral is: 



W(l) = (0(ri)0(r 2 )0(r 3 )r 13 |T fil (l) +^,(1) + f Jj3 (l)|0'(r 1 )0 / (r 2 )0 , (r3)ri 2 ) (38) 

The radial part can be separated into three contributions. The resulting integrals can be 
straightforward evaluated using the operators Eq. (14-16) and the definition of the two- 
electron integrals Eq. (A. 3): 
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j /,\ _ ( n l + n l) t/ at 9 /vr at.,, ,, i \ "ii ,m'i ,mi ,m 2 ,m3 ,m 3 /o Q \ 

Jxe,Ri(iJ g wi,w2,w 3 ,i,i; Ji)J , ) , 2) , i) , 3) ,, (dyj 

+ {n l u l + u l )J{N 1 -l,N 2 ,N 3 ;u l ,u 2 ,u 3 ;l, ±) w l w ? 

-^ J (^i^2,^ 3 ;u;i,a; 2 ,W3;l,l) Zl/i ^ w J 3 , (40) 



7" m — ( n l + ^ 7V at at at-.,, ,, ... -i 1 \miy i ,ui2,m> 3 ,m3 

JxE,R 2 (lJ ^ J(7Vi,7V 2 ,7V 3 ,a;i,a;2,W3,-l,i;, i) ,, ) , a) j, ) , S) ,, 

+ y J l 7V i + i ' 7V 2,-'V3,Wi,a;2,W3,-l,i;, li i i)j2i , i) ,3 )J / 

(J- - n l) t/ at 9 at , 9 at.,. „ ,,. i 1 ^i,m' 1 ,m2,^,m3,m 3 
J^Vi - Z,l\ 2 + l,lSz,Ui,U2,Uz, -1, i ;{ 1 ,i' 11 i a ,^,i 3) iJ 



- fj( Nl - 1, JV 2 + 2, tf s ; u*, u*, c 3 ; -1, l)D5S5' m3 '^ ( 41 ) 
Note the symmetry of these integrals: 

J(Aq, N 2 , N 3 - u u u 2 , u 3 ; -1, l)T^JST Z,m ' 3 = J ^ N *> N * ^ "*> l )2^T 2,m ' 2 

(42) 

The last term: 

W 3 (l) = W^Nmws; l,l)SlS7r m3 '^- (43) 
Now let us evaluate the angular kinetic energy. The integral containing the square of the 
angular momentum operator can be easily calculated by means of its eigenvalue equation: 



Ike^l = ^(r 1 )0(r 2 )0(r 3 )r 1 3|^|0 / (r 1 )0 / (r 2 )0 / (r 3 )r 12 ^ 

= + 1) J( Nl - 2, N 2 , N 3 ; ^ c 3 ; 1, ( 44 ) 

where l[ is the quantum number of (f>(ri). 

The evaluation of the other angular kinetic energy contributions is more involving. The 
following integrals should be evaluated: 

lKE,e 1: i = -{ 0(ri)0(r 2 )0(r 3 )ri 3 |— —^j^T, — |0(ri)0(r 2 )0(r 3 )r 12 ) , (45) 
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iKEjhfl = ~\ ^(ri)0(r 2 )0(r 3 )r 13 l ^ ^- ^ cot Oi^r |0(ri)0(r 2 )0(r 3 )r 12 ^ , 

(46) 

^*2 sm ^^2 (-^^ 

^,0 1; 3 = - (0(ri)0(r 2 )0(r 3 )ri3|— — - sin (0i - <h)^r^ — |0(ri)0(r 2 )0(r 3 )ri 2 ) . 

(47) 



3.1 Evaluation of Ike,0i 

The integral 



I ^*2 COS ^2 \ 

( 0(ri)0(r 2 )0(r 3 )ri 3 | — I0(ri)0(r 2 )0(r 3 )r 12 ) = 

\ r\Tyi sin 9i aviary / 

L 3=\h-l 3 \ 



X 



Jo Jo 

Y^(e 3 ,<f> 3 )<m 3 



x 

7T p2lT 



X 




^0 



cos^r™* 2 (^ 2 ,0 2 )F™ 2 (^ 2 ,0 2 )^ 2 



7T /»27r 1 ^ 



X X X ^ ^'^VW.W^ (48) 
with c 3 l 3 = C Ls (l 3 , m' 3 ; / 3 , m 3 ). dQiis the volume element sin(6 l i)(i^i(i0i. Let us perform first 
the angular integration. The derivatives of a spherical harmonic Y™ 1 (0 1 , <pi) with respect to 
the polar angle 6>i [? , 5.7.] is: 

dY?HB 'i,0i) 

where: 



f la e-^Y i r^\9 1 , <P,) - f lb e^Y [ r-\9 1 , 00, (49) 



/i« = ^'1 + m[ + 1)(Z1 - ml)] 1 / 2 , f lb = - mi + l)(Jj + m',)] 1 / 2 (50) 

multiplying by the complex conjugate Y™ 1 * (9i, <fii) and linearizing the product using Eq. 
(A. 4) we obtain: 
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E [2^ + 1 ] 



1/2 



(4tt) 1 /2 

x {fiaCiaY^ 1 ^,^-^ - f lb c lb Y^- 1 (e 1 , ( p 1 )e^} (51) 



with 



ci a = C Ll (Zi,mi + l;Zi,mi), Ci 6 = C Ll (Zi,mi - l;Zi,rai). 
On the other side using the cosine recursion relation: 



(52) 



cos 6 2 Y^ (9 2 ,<h) = f2aY;p 1 (9 2 ,h) + f2bY l 'ri 1 (6 2 , 2 ), 
for Z 2 — 1 > 0, otherways the first terms vanishes. The factors are: 



2a 



(l' 2 + m 2 )(l 2 -m' 2 ) 
[ (2l' 2 + l)(2l' 2 -l) 



1/2 



f2b — 



(l' 2 + m' 2 + l)(l' 2 -m' 2 + l) 
(2/ 2 + l)(2Z 2 + 3) 



1/2 



multiplying by Y t ™ 2 * (Q 2 , $2) and linearizing the products: 



1^*(02,&) COS 02^(02,02) = 

? — - J] [2L 2 + I] 1 / 2 faCaYgfa, fc) 
^ ' L 2 =\l'o-l-h\ 



+ 



1 



(53) 



(54) 



(47r)V2 Z 

V 7 L 2 =|« 2 +l-Z 2 | 



E [24 + l] 1/2 /2 b c 2b F i 5(0 2 ,0 2 ) (55) 



with 



c 2a = C L2 (Z 2 - l,m' 2 ;l 2 ,m 2 ), 

The rotation of Y™ 2 (d 2 ,4> 2 ), see Eq. 
47rY^ 2 (0i,0i). Altogether (no factor 4pi): 



c 2b = C L 2(Z 2 + l,m 2 ;/ 2 ,m 2 ). 
(A. 7), and integration over 12 leads to 
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1/2 



E E [(2^i + l)(2L 3 + l)] 1/2 

Li=|/i-Ji| ^3 = 1^— « 3 | 

V 17 [L 2 =\l' 2 -l-h\ 

{haC la f2aC2aY^ + \6 l , ^)Y^ (9 U ^)e^ 

-fi b c lb f 2 aC2aY^- 1 (9 1 , fa)Y%\O u fae**) 

l' 2 +l+h 

+ E [ 2L 2 + X ] 1/2 (/laCla/ 2fe C 2 ^ 1 + 1 (^, h)Y^(9 1 , 0^"^ 

L' 2 =\l' 2 +l~l 2 \ 

-fi b ci b f 2b C2 b Y^- 1 (9 1 , fc)Y$(O u 0!)e^) } (56) 
Now the spherical harmonics with the same arguments are linearized again, see Eq.(A.4). 

c 3L = C L3 (l' 3 ,m' 3 ;h,m 3 ). 
Note we need the complex conjugate in Eq. (A. 10). (easier would be (-1) "M2 ) 

(An)^(-l)^ E E [(2^ + l)(2Z 3 + l)] 1/2 c 3i3 ("-^) 

{l' 2 -l+h L!+L 2 
E E (-l) M2 [(2£ 2 + l)(2L + l)p 

L 2 =\l' 2 -l-l 2 \ L=\L 2 -Lx\ 

(flaC la f2aC2aC, a Y L f(9 1 , <h)Y* +M > +1 (9 U ^e"^ 
l' 2 +l+h Lx+L 2 

+ E E (-l)^[(2^ + l)(2L' + l)]^ 

L' 2 = \V 2 + l-h\ L'=\L' 2 -Li\ 

{flaClaf2 b C2 b C4cY^(9 1 , < f> 1 )Y^ +M ^ 1 (9 U ^)e~^ 

-fi b ci b f2 b C2 b C4 d Y^(9 l7 1 )V L Ml+M2 ~ 1 (^i, 0i)e^) } (57) 
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with 



x 



c 4a = C L {L X , M 1 + 1; L 2 , -M 2 ), c ib = C L (L U M x - 1; L 2 , -M 2 ) 

c 4c = C L \L 1 ,M 1 + l ] L' 2 ,-M 2 ), c 4d = C L '(L 1 ,M 1 -l;L' 2 ,-M 2 ) (58) 

(47r)1/2( _ 1) M 2+ M 3 £ £ [(2L 1 + 1)(2L3 + 1)] 1/2 C3,3 --^ 

L 1 =|J' 1 -J 1 | L 3 =\l' 3 -l 3 \ ^ X ' 

1' 2 — l+Ja Li+L 2 

E E E (-l) M2+M3 [(2L 2 + l)(2L + l)(2L 4 + l)] 1 /2 

L 2 =|Z^-l-Z 2 | L=|L 2 -Li| L 4 =|L-i 3 | 

(/uCl./ a aP a .C te C BB yJ 1+Ma+ ^ +1 (fl 1 ,0 1 ) C -^ 

-/l6Cl 6 / 2 aC2aC4 6 C 56 l^ 1+M2+M3 - 1 (0 1 ,(/» 1 )e^) 
Z 2 +l+/ 2 Li+L 2 ^+^4 

+ E EE (-l)^ +M3 [(24 + l)(2^' + l)(2L' 4 + l)] 1 / 2 

L 2 =|/ 2 +l-Z 2 | L' = |L' 2 -Li| I4=|L'-2, 3 | 

(/iaC la / 2b c 2b c 4c c 5c F^ +M2+M3+1 (^ 1 ,0 1 )e- i ^ 

-/uC! J 2b c 2b c 4d c 5(i r* /l+M2+M3 - 1 (^i, 0i)e^) } (59) 

needed are 

c 5a = C L *(L,M 1 + M 2 + 1;L 3 ,-M 3 ), c 5b = C U (L, My + M 2 - 1; L 3 , -M 3 ) 

c 5c = C L H^,M 1 + M 2 + 1;L3,-M 3 ), c 5d = C L '*(L', M 1 + M 2 - 1; L 3 , -M 3 ) (60) 

Finally, the recursion relation containing the sine function is applied [14, C7]: 

"(2Z + 1)^ 1/2 



sin 6^ 2m 



(2/ + 3) 



[[(/ - m + 2)(Z - m+ l)] 1 ^ 1 ^™- 1 ^, 00 



+ ((/ + m + 1)(Z + m + 2)) 1 /V^ 1 F™ 1 +1 (# 1 , 00] , (61) 

The relation is valid for m^O. This condition is fulfill in our cases M = Mi + M 2 + M 3 + 1 
and M' = M X + M 2 + M 3 - 1. Multiplying l/sin0i and Eq. (40) we have: 
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li+h l' 3 +k 

I K E, ei (l) = W 1/2 (-l) M2+M3 E E [(2^ + l)(2L 3 + l)] 1/2 c 3 ,3 

Li = |«i-/i| L 3 = l«^-i3l 

l'2-l+h L1+L2 T/or , i\/nr , 1 M 1/2 



with 



E E E (-i) M2+M3 (2^+i) 

/, 2 =|^-l-Z 2 | L=|L 2 -Li| L 4 =|L-L 3 | 

(*7T / i 27T 



(2L 2 + 1)(2L 4 + 1) 



(2L + 3) 

^ 1 1 +M2+M3 (^l,0l)^l 



F Mi+M 2+ M 3+ 2^ i ^ i)e -2^ 1 ^ i 

n27T 
r M 1 +Af 2+ M 3 -2^ i)0i)e 2^^ i 



x - 1, N 2 + l,^;^,^,^; -1,1; L 2 ,L 3 ) 

'(2L 2 + l)(2L / 4 + l) 



/ 2 +l+«2 Li+L 2 r/nr , lX/rir/ , n \-il/2 



+ E E (-i)^ +M3 (2£'+i) 



L 2 = |« 2 + l-i 2 | L'=|L 2 -Li| 



(2L + 3) 

/»7T pzit 
flaClaf2 b C2bhcC, c C 5c / / 

Jo •/ 



2 — I'jTi" '2| ^ —|-^ 2 

7T />27r 



'0 JO 

r7T /*2n 



F Mi+M 2+ M3 + 2 ( ^ )0i)e -2^ 1 ^ i 

n27T 
F ^l+M 2 +M3-2^ i)0i)e 2^ 1 ^ i 

n2"7T \ 
Y%£ M * +M >(O u <j>M dQ, 

xJ(N, -1,N 2 + 1, iV 3 ; Wl , w 2 , w 3 ; -1,1; £ 2 , £3)} (62) 



/3a=^p-M + 2)(L-M + l)]^ 



/36= 21? [(jL + M + 2)(jL + M + 1)]1/2 ' 
h c = ^\{L-M' + 2){L-M' + l)f\ 

h d =^[(L + M' + 2)(L + M' + l)}^ (63) 



2M 

The first kind of integrals over spherical harmonics are: 
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7T / 1 2-7T 



JO 



M1+M2+M3, 



Lj 0i)tK2i = S(L + 1, 0)(J(Mi + M 2 + M 3 , 0) = 



(64) 



These integrals are zero because L > 0. The following integrals are evaluated integrating 
over the associate Legendre functions : 



A Y1 (L, M + n, -n) 



7T p2tt 



jo 

tt 1 / 2 (2L + 1) 1 / 2 



(L — M — n)\ 



1/2 



P^ /+n (cos0i)sin0id0i (65) 



(L + M + n)!_ 

using the algorithm developed by Wong [18], in which the overlap integral over Associate 
Legendre functions is: 



Pi max Pi max 



(cos0)P L A f(cos0)sin(0)d0 = £ E^U^m, 

Pl=0 P2=0 

T (|(Li + L 2 - M x - M 2 - 2p! - 2p 2 + 1)) T (i(Mi + M 2 + 2 Pl + 2p 2 + 2)) 



(66) 



r(i(Lx + L 2 + 3)) 

T are Gamma functions. p max = [(/ — m)/2] the integral part of (I — m)/2. This formula is 
valid for < M\ < L\ and < M 2 < L 2 . For M\ < L\ or M\ < L\ the integral is zero, see 
Ref. [18J. If M is negative the formula is used: 



p-M, 



COS 0i ) 



^ {L - M) \ P f[co^). 



L " v (L + M)\ 

For even L\ + L 2 — Mi — M 2 the integral vanish. The coefficients are: 



a, 



■1)p(L + M)\ 



lm 2 m + 2 p(M + p )\p\(L- M -2p)\ 
The radial three-electron integral fulfill the symmetry : 



(67) 



(68) 



J(JVi - 1, N 2 + 1, N 3 ; wi, w 2 , w 3 ; -1, 1; La, £3) = 

J(iY\ - 1, N 3 , N 2 + 1; wj, w 3 , w 2 ; 1, -1; L 3 , L 2 ) 



Finally the programmable expression is: 
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l[+h l' 3 +l 3 ( l' 2 -l+h L x +L 2 L+L 4 

w(i) = ( fc )'/ 2 E E E EE 

Li=|«i-Zi| Ls=|^-/ 3 | [l 2 =|^-1-/ 2 | L=|L 2 -ii| i 4 =|i-i 3 | 

(~ 1 ) M2+M3 (|^p K 2L i + !)( 2L 3 + 1)(2^2 + 1)(2L 4 + I)] 1 / 2 

UlaClaf2aC2afzbCzL A C Aa C ba A Y l{L + 1, M + 1, -2) 
~ hbClbf2aC2ahcCzL 3 C ih Cr >h A Y l{L + 1, M' — 1,2)) 

x J(Aq - l,^3,^2 + l;wi,u;3,a;2;l,-l;^3,L 2 ) 

2 2 +l+« 2 Li+L 2 L'+L 3 / or . -| \ 

, \" / -|\M^+M 3 

+ 2^ ^ ^ (2L + 3) 1 /2 

L 2 = |« 2 +l-« 2 | L'=|L 2 -Li|L^=|L'-L 3 | V ' 

[(2Lx + 1)(2L 3 + 1)(2L' 2 + 1)(2L; + 1)] 1/2 

(fla c laf2bC2bf3bC3L 3 C4 :C C 5c AYl (L' + l,M + l,-2) 
-flbClbf2bC2bhcC2,L 3 C Ad C bd A Y l{L' + 1, M' — 1, 2)) 

x J(Aq - 1, iV 3 , iV 2 + 1; Wl , W3 , w 2 ; 1,-1; L 3 , L' 2 )} (69) 



3.2 Evaluation of 



(92 



For the evaluation of I re, 92, the same steps than for Ike,$i Eqs. (31-34) are needed. 
Afterwards the product of spherical harmonics of electron 2 should be linearized. Then a 
rotation like in Eq. (A. 7) is performed and the complex conjugate of Y^ 1_1 (#i,0i) used. 
The new products of Y^ Ml+1 *{9 u fa)Y% 2 fa, fa) lead to: 



Y™s*(e 2 , 2 )E ; r% 2 , &)) <wr%, 2 )) 



3V m i/fl ^\\ i 2 +(2 * 3 +*3 

%/ U(4,)n-i)^". e E E 



Li+L 2 i+i 3 

^ ^ [(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L + 1)(2L' + 1)] 1/2 c 2L2 c 3L , 

L=\L 1 -L 2 \ L'=\L-L Z \ 

x { h a c la c 2a Y^+ M ^ + \e u fae-* 

-fibC^bY^^-'^^e*} (70) 
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where 



c 2 l 2 = C L2 (l' 2 ,m' 2 ;l 2 ,m 2 ), c 3L . A = C L3 (l' 3 ,m' 3 ;l 3 ,m 3 ) 
c 2a = C L (L 1} M 1 + 1; L 2 , -M 2 ), c 2b = C L (L 1} M 1 - 1; L 2 , -M 2 ) (71) 
c 3a = ^'(Li, M x + 1; L 3 , -M 3 ), c 3b = C L '(L 1 , M l - 1; L 3 , -M 3 ) (72) 

Now we apply the recursion relation including the cot6> function [? , 5.7.]: 

-cot^r%,0i) 



which leads to: 



-[{I + m + 1)(Z - m^V^E™ 1 ^, 00 



2m; 1 

+ 2^7 [(/ " m + + m )] ll2el<t>lY r~ 1 ^ ( 73 ) 



^i+Zi Z^,+Z 2 Z 3 +J 3 Li+L 2 £+£3 

w, 2 = v^(-i)* + «- E E E E E 

L 1 = \l[-h\ L 2 =\l' 2 -h\ L 3 =\l' 3 -l3\L=\L 1 -L 2 \L'=\L-L 3 \ 

[(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L + 1)(2L' + 1)] 1/2 c 2L2 c 2L3 

r>7T /-27T 

X 



nZ7T 
F M 1+ M 2+ M 3+ 2 ( ^ )0i)e - 2l , 1 ^ i 

n27T 
y M 1+ M 2+ M3 ( ^ 0i)dni 

n2"7T 
r M 1+ M 2+ M 3( ^ ?0i) ^ i 

Y M, + M 2+ M. A -2 {6iAi y t4>ld ^ 

x [J(iVi - 2, N 2 , N 3 ; Ui, u 2 , 003, 1, 1; L 2 , L 3 ) 
- J(N 1 ,N 2 ; N 3 ; ooi,oo 2 , oo 3 ; -1, 1; L 2 , L 3 ) 

- J(JVx - 2, iV 2 + 2, iV 3 ; Wl , W2 , w 3 ; -1, 1; L 2 , L 3 )] , (74) 

with M = M 1 + M 2 + M 3 + 1 and M' = M 1 + M 2 + M 3 — 1. Integrating: 
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ha=^[{L' + M + l){L'-M)Yl\ 

f2 b =^W-M+l){L' + M)Y/\ 
f2c = ^[(L' + M' + l)(L'-M')]^, 

f2 d = ^[(L' + M' + l)(L' + M')] 1 / 2 (75) 

Using Eqs. (46-49), and the symmetry property of the radial three-electron integrals, the 
final programmable expression is: 

l[+h l' 2 +h l' 3 +h L ± +L 2 L+L 3 

w^ac-d""* E E E E E 

Li=|/i-Ji| L 2 = \l' 2 -l 2 \ L 3 =\l' 3 -l 3 \ L=\L X -L 2 \ L'=\L-L 3 \ 

[(2L X + 1)(2L 2 + 1)(2L 3 + 1)(2L + 1)(2L' + 1)] 1/2 c 2L2 c 2L3 

X {flaClaf2aC2aC3aA Y l( y L', M+1,-2) 
+ ^(L',0)^(Mi+Af2+M 3 ,0)(/laCla/2bC2aC3 a — flbC\bf2cC2bCzb) 
~ flbClbf2dC2bC3bA Y l(L', M-1,2)} 

x [J(^ - 2, 7V 2 , 7V 3 ; Wl , W2 , w 3 ; 1, 1; L 2 , L 3 ) 
- J(N 1 ,N 2 ; N 3 ;uj!,uj2,uj 3 ] -1, 1; L 2 ,L 3 ) 

- J(7Vi - 2, iV 2 + 2, iV 3 ; Wl , W2 , oo 3 ; -1, 1; L 2 , L 3 )] , (76) 

3.3 Evaluation of Ike,4& 

The derivative of the spherical harmonic with m\ = with respect to <pi vanishes: 



Wi,3 = 0. (77) 

In case of mi 7^ 0: 

S 2^=imlHM). (78) 

In this section we evaluate Eq. (30). The function sin(0i — 2 ) can be written in exponential 
form: 
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8111(0! - 2 ) = i (e^e-* 1 * - e-^e^) , (79) 



and the following functions can be written as spherical harmonics: 

sinfce"^ = y^yf 1 ^,^), sin^e^ 2 = -^Y 1 1 (e 2 ,h)- 
The products of spherical harmonics can be linearized: 



(80) 




iri^Yr 1 ^) = yf E ( 2L2 + 1) 172 ^ 1 "" 12 ^,^) 

L 2 =|/2-l| 

^(^,02)^(^2,02) = £ (2L 2 + l) 1 / 2 c 2b Yl;^(e 2 ,<j )2 ) (81) 




L2 = |«2-l| 

with 

c 2a = C L2 (l,-l;Z 2 ,m 2 ), c 2b = C L2 (l,l;l 2 ,m 2 ). (82) 
Using the above derived expression, the product of angular functions is: 



sin On dY^iOiAi) 
-Y l ^(0 1 , 0O>T(^ ^)— ^ sin (0, - 2 )^^ Y^(0 2 , 2 ) 

1 mi /2 ^ 3 h+1 



( 3+'3 (2 + 1 

3 E E [(2^2 + 1)(2L 3 + 1)] 1/2 cs^C 3 ^, 0a) 



(4tt)V2 2 

V ' L 3 =\l' 3 -l 3 \ L 2 =\h-1\ 



pt<Pi 

-l-m 2 /n ± \v m l(Q A \ ~" m ^ 



with C3L3 = C L3 (/ 3 , m' 3 ; /3, 7713). Using the recursion relation containing the sinus function 
Eq. (42) over Y™ 1 {61, (pi) (this relation can be always used directly over Y™ 1 (61, 0i) because 
by definition m' 1 7^ 0, for m' 1 = 0, the derivative of a spherical harmonic is zero, see Eq. 
(xx)) and the factors defined as: 
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f la = [(l[-m' 1 + 2)(l[-m' l + l)] 1 / 2 , 
Eq. (63) can be written: 



/ lb = [(/ / 1 + m / 1 + 2)(/ / 1 + m / 1 + l)] 1 / 2 



(-l) m 2 



1 1 2 



(4tt)V2 4 V 3 



2l[ + 1 



1/2 l 3+ l 3 l 2 + l 

E E [(2L2 + 1)(2L S + 1)] 

L 3 =|^-/ 3 | L 2 = |/ 2 -l| 



1/2 



_2l[ + 3_ 

X C3L 3 ^(^,03) {c 2 a^(^2,02)y L - 2 1 " m2 (^,02) 

+ C 2 5^(^2,02)V L 1 2 - m2 (^,02) 
7la^(^l, M^r^l, + /l^r ^' Me'^^FWlM) } 

Combining F^(^,0 2 )F L2 1 - m2 (^,0 2 ): 







—i) m ' 2 e ^l 1 /! 172 ^^ 1 ^-^,) 



l 2 =|l 2 -; 2 | 



(4tt)V2 



with 



c 2 L'a = C L ' 2 (L 2 , -1 - m 2 ; Z 2 , -m 2 ), c 2L 'b = C L ' 2 (L 2 , 1 - m 2 ; Z 2 , -m' 2 ) 
Eq. (65) can be rewritten as: 



,> 1 1 /2 



4tt4 V 3 



07 , . -,-,1/2 *3+*3 2 2 +l L 2 +l' 2 l'i+1+h 

Ll-y + 1 V — * \ — * \ — * \ — * 

21' -f 3 ^— ^— -/ ^-^ 

1 J L 3 =|^-« 3 |L 2 =|Z 2 -l|L 2 =|L 2 -i 2 |Li=|i' 1 +l-Zi| 

x [(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L 2 + 1)] 1/2 

X {c 2a C 2L , a C3 L 3r L f S (^3,03)^ 2 ~ 1 (^2,02) 

x (ZiaC^-^^^Oe 2 ^ +/ 16 c 16 y^ +1 (^ 1 ,0 1 )) 

+ C 2fc C 2L , fe C3L3^ 3 (^3, 3 )^ 2+1 (^, 

{fiaCiaYg-^, 00 + frttoY£+\B u 0i)e- 2 ^) } , 
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with 



c la = C Ll (l[ + 1,771; - l;Zi, mi), c lb = C Ll (l[ + l,m; + l;Zi,mi). 
After a rotation of the spherical harmonics of electron 2, see Eq. (A. 7): 



(-I)' 



4tt4 V 3 



2l[ + 1 
2/; + 3 



-| 1/2 J 3+'3 «2 + l 



E E E E 

L 3 =|^-«s| L 2 =|«2-l| L 2 =|L 2 -i 2 | Li=|ii+1-Ji| 



X [(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L' 2 + 1)] V2 C 3 L s Y£ S (03, fo) {c2 a C 2 L>a 

fiaC la J Q J Y^- 1 (9 1 , ( f )l )Y^- 1 (9 1:( f )l )e 2 ^dn 1 

r L ^ +1 (ei,0i)y L Y 2 " 1 (^i^i)^i) 

flaCla J Jo y ^ 1 " 1 ( 01 '^ 1 )^ a+1 (^i^i)^i 

n27T 
E L ^ +1 (^i,0i)E L Y 2+1 (^i,0i) e -^Mni 

rotating Y^(6 3 ,(f> 3 ) we obtain 47rE L f (^, 0i): 



*i,S = (-1) 



j 1 /2 



4 V 3 



2/; + i 



2l[ + 3 



1/2 *3+'3 Z 2 + l 



L 2 +* 2 



ii+l+Jl 



E E E E 

L3 = \l^-h\ L 2 =\l 2 -1\ L' 2 =\L 2 -l' 2 \ Li=|ii+l-h| 

x [(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L 2 + 1)] 1/2 {c 2a c 2L , a C3L 3 

x(/laCi a ^ jf r^^i^i)^- 1 ^,^)^ 2 " 1 ^^!)^ 1 ^! 

r»7T /»27T 



nZ7T > 

(/*7T /*27T 



7T /-27T 



+ /lbCifi 




JO 



taking the complex conjugate and using the product theorem: 
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l KE,<t>i,3 



8 V 3tt 



07 , . t -i 1/2 '3+ ; 3 h+l L 2 +V 2 l' L +l+h 

1 E E E E 



2/; + 3 



L 3 = |Z! i -«3| L 2 = |/ 2 -l| L' 2 = |L 2 -«^| ii=|«i+l-ii| 
L 3 +L' 2 

[(2^i + 1)(2L 2 + 1)(2L 3 + 1)(2L' 2 + 1)(2L' 3 + 1)] 1/2 {c 2a c 2L , a c 3L3 

Z^ = |L 3 -L^| 

fiac la c 3a J o J Y^- 1 (e 1 , ( f )1 )Y^ +M ^ 1 (e 1 , ( p 1 )e^dn 1 



F L ^ +1 (^ 1 ,0 1 )^ 2+M3 " 1 (^i ! 0i)^i 

/iac laC3a y o y o y^ i-i (#i, 0i)y^ 2+j ^ 3+i (^i, 

n27T 
r^ +1 (^ 1 ,0 1 )^ 2+M3+1 (^^i)^ 1 ^i 
3 



with 



(92) 



c 3a = C L 3(4, M 2 + 1; L 3 , -M 3 ), c 3b = C L '^L' 2 , M 2 - 1; L 3 , -M 3 ). (93) 
And defining the new auxiliary angular integral: 



P^ 1+n (cos 6)P^ +n (cos 9) sin fldfl (94) 

This integral of the product of Legendre polynomials is solved according the Eq. (47) from 
Wong, finally the integral is: 
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1/2 l 3+ l 3 h+l L 2+ l 2 li+l+h 

E E E E 

L 3 =\l' 3 -l 3 \ L 2 = \h-l\ L' 2 =\L 2 -l' 2 \ ii=|«i+l-ii| 

^ [(2Li + 1)(2L 2 + 1)(2L 3 + 1)(2L' 2 + 1)(2L' 3 + 1)] 1/2 {c 2a c 2L , a c 3L3 

x (/i„ci a c 3a A y2 (L 1 , Mi - 1, L 3 , M 2 + M 3 - 1, 2) 
+/i 6 c 16 c 3b (-l) Ml+1 <5(L 1 , 4)(J(M! + 1, M 2 + M 3 - 1)) 

+ C 2b C 2L , b C 3L:i {fiaCiaCsai-l)^-^^, L' z )5{M l - 1, M 2 + M 3 + 1) 

+ fibCibCsbAyziL!, Mi + 1, L' 2 , M 2 + 1, -2))} (95) 

We have developed a computer code and computed the radial integral equations Eq. 
(24-26) and angular integral contributions to the kinetic energy Eqs. (27,50,56,72). The 
computer code has quadruple precision. In Table 1 we show some selected values of two- 
electron kinetic energy integrals calculated with high precision. The integrals of Table 1 
have been checked with the ones calculated using the Kolos and Roothaan transformation, 
agreeing to more than 30 decimal digits. Finally, in Table 2 several kinetic energy integrals 
have been decomposed in all their components. 
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Appendix A: conexion to KR method. Why the KR method cannot be extended 
to three-electron integrals 



The derivative of r\ 2 is: 



lT = -V 2 r 12 = Vi 2 r 12 = ur 12 — 

r i2 



Short proof: 



(A.l) 



V 1 =V 1 = — i +t—J +T^k 
ox i oyi az\ 

the interelectronic distance in Cartessian coordinates: 



(A.2) 



ri2 = V {xi - x 2 ) 2 + (yi - y 2 ) 2 + (zi - z 2 ) 2 



(A3) 



dx\ dx\ 



dx\ 
1 

ri2 



[xi - x 2 )~t + (y 1 - y 2 )~i + {zi - z 2 )~t 



T\ 2 



ri2 
ri2 



(A.4) 



The scalar product of a vector with itself: 



ri2 _ rn 
r\ 2 r 12 



\ r 12\ 



12 
12 



r 2 - r 2 - (A.5) 

'12 ' 1 

In the case of the two-electron kinetic energy integrals we can make use of the property 
Eq. (10) and to obtain the expressions in terms of radial coordinates r 12 : 



dri 



r\2 r\2 



+ J r\r' [Vi0Kri)] [VMr^r^'dn 
J r 12 



+ u' J [^VMr^r^'-'—dn (A.6) 
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The three-electron kinetic energy integral: 



<^r 2 )0' 2 (r 2 )dr 2 \ 0*(r 3 )0 3 (r 3 )dr 3 [ [<PUr 1 )r 12 V* 1 ][V 1 r 13 ^ 1 (r 1 )]dr 1 (A.7) 



using the turn over rule of the antihermitian operator nabla: 

V* = -V (A.8) 

I = J [0t( ri )r 12 Vt] [Vir^fa)] dr x (A.9) 

the derivative of a product: 

Vl0;(ri)ri 2 = r 12 Vt0K r i) + 0!( r i)Vlr 12 

Viri 3 0i(ri) = ngVi^UrO + / 1 (r 1 )Vir 13 (A.10) 



[0!(ri)r 12 Vt] [Vir^UrOJdn^ 

ri 2 ri 3 [V^(ri)] [V^rO] dn + / ^r^ifa) [V^r 12 ] [Vir 13 ] dn 
+ y [0Kri)V 1 / 1 (r 1 )][r 13 V 1 r 12 ]dr 1 + | [^(rOV^M [r 12 Vir 13 ] dr, (A.ll) 
Let us evaluate first the integral f 0*(r 1 )0' 1 (r 1 ) [V*r 12 ] [Vir 13 ] dri: 

/ ^(rO^'^n) [V;r 12 ] [Vir 13 ] dn = / #( ri )&(n) — • — ^ (A.12) 

J J ?"l2 ?"l3 

the scalar product of the two vectors rj^ and leads to: 



ri2 • ri3 



(xi - x 2 )~t + (y 1 - y 2 )~t + (zi - z 2 )~t 



x 



(xi - x 3 )T + (yi - y 3 )Y + (zi - 2 3 )Y 



(xi - x 2 )(xi - x 3 ) + (yi - y 2 )(yi - y 3 ) + {z x - z 2 )(z x - z 3 ) (A. 13) 
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For simplicity, let us consider only the coordinates x: 



(x l - x 2 )(x 1 - x 3 ) —x\- xix 2 - + x 2 x 3 (A. 14) 

together with the coordinates y and z, it can be transformed to: 



ri2 ■ rla = ^ {r\ 2 + r{ z - r% 3 ) (A.15) 



so the integral is transformed to: 



I MnWn) [v;r 12 ] [V iri3 ] dn = J fi(nWri) r ^ ■ — dn = 



+ 2 



ri2 r 13 
- / 0t(ri)0i( r i)^ rfr i 

1 J r 13 

\ [ MnWn^dn-l [ ^^-^-dn (A.16) 

^ J T\ 2 I J T\ 2 T 13 

Note that a new integral type has appeared, namely: 

/ -^-\ (A.17) 
\r\ir\ 3 l 

Comparing with the case of the two-electron kinetic energy integrals, where a simple overlap 
integral was obtained, in the case of the three-electron integrals, the derivatives of the 
interelectronic coordinates produce an scalar product of two vectors, which is transformed 
in a more complex relation including interelectronic distances. Problem is the new integral. 
Now we study the next integrals: 

j [ftOrOVi^rO] [r 13 Vir 12 ] , j W^MPi^i)} [^V^] dn (A.18) 
making use of the rule of a product: 

Vi [r 12 r 13 ] = r 12 Vir 13 + r 13 Vir 12 (A. 19) 

we add and substract, for instance, the integral: 

J [^(rOV^i^i)] [ri 2 Vir 13 ] dn (A.20) 

combining 
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^(rOVi^rO] Vi [r 13 r 12 ] d n = J [^(r^VM^)] [r 13 Vir 12 ] dn 

+ J [(t>* 1 (r 1 )V 1 (l ) ' 1 (r 1 )}[r 12 W 1 r 13 }dr 1 (A.21) 



considering 



J [0J(ri)Vi0i(ri)] Vi [r 13 r 12 ] dn = -J r 13 r 12 Vx ^(rOVi^rO) dn (A.22) 
according to the Green's identity: 

- V (^V0) = -^V 2 - (W>) (V0) (A.23) 

combining: 



ri 2 r 13 [V^KrO] [Vi^rO] dn - J r 13 r 12 Vi (#(ri)Vi#(n)) ^ = 

- y ^^^(rOV^UrOdn (A.24) 

AUtogether we have: 

- y [^(rOnaVI] [Vin^UrOJdn^-ynW^rOV^UrOdn 

* J r±2r 13 



2 

+ 



y [0'i(ri)V^t(ri)] [r 12 Vir 13 ] dn 

- y [^(rOVi^^rOltnaVirxaldn (A.25) 

From this point we cannot proceed as in the two-electron Kolos and Roothaan transfor- 
mation, because it is not possible to put together the interelectronic distances into a common 
derivative, with the exception of the steps done already Eqs. (34-39). 

Other way to solve the last two integrals is to evaluate them individually. For instance, 
let's take: 
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J [^(rOV^Kri)] [rigVina] tZn = J {V\(j)\{r x )) (Vir 12 ) dn (A.26) 

The derivative of an orbital, see Ref. [26], Appendix D: 



ih"(»i,0i) (A.27) 



V *( r '»=(^ t + ^ + ^ t )W'» (A ' 28 > 
Let's consider the derivative with respect the coordinate x and use the chain rule of differ- 
entiation: 

<90i(ri) _ dri dcpiiji) t d6 1 <90i(ri) + dfa dfajrx) ^ 



dx\ dx\ dr\ dx\ d6\ dx\ 
The terms with coordinate r\. 



(V^Kn)) (Vir u ) 

dri-} dri — > dr 1 ->\ ( <9r 12 — >• dr 12 — > dr i2 -^\ dcbAr-C, 
1 + ~r\ J + ^— « -t; — t + -t; — J + -t; — « 



ctei (9zi / \ dxi dyi dz\ ) dri 

dri dri2 dri dri 2 dri dri 2 \ d(pi(ri) 



dxi dxi dxi dx\ dx x dxi J dri 
xi (xi - x 2 ) yi (yi - y 2 ) z x [z x - z 2 )\ <9</>i(ri) 



r\ ri2 ri ri2 n r i2 / dri 

r i2 + r i — r i <90i(i% 



2rir 2 9ri 



(A.30) 



The factor: 



ri2 + r\-r 2 2 
2r x r 2 

appears in Ref. [26] Appendix C and Ref. pL4j Eq. (55). 



(A.31) 



r 12 0' 1 (r 1 )(V^ 1 (r 1 ))V 1 r 13 dr 1 = /V 12 ^ri)^^ — -dn + ■ ■ ■ = 

J dri r 13 n 

/. Wl(r0 ^ri±r^i, ri + .., A , 2) 
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The term corresponding to the variable 9\. 



(V;#(rO) (Viru) 



d6i-$- dOi-t dOi -A / 0r X2 -> <9ri 2 -^ <9ri 2 ->\ <90i( r i' 

2 + 7T~ J + 7^— K "Tj 1 + "7S J + -Ti K 



9xi 9yi / \ 9xi ch/i (9^i / <96*i 

<9r 12 ^ 0r 12 ^ 09 1 0r 12 \ 0<fii(ri) 



0x\ 0x\ Oyi 0x\ 0z\ 0x\ J 09\ 
This expression can be transformed into polar and interelectronic coordinates, as shown in 



Ref. p2], with: 



cos^) = , 2 Z \ 2 (A.33) 



09 x 


Z1X1 




Ox 1 


rf sin 9\ 




09 x 


Z1V1 




Oyi 


rf sin 6*i 




09t 




1 


0z\ 


rf sin^! 


ri sin#i 


as Ref. 


[15]: 





(A.34) 



09i0r 12 09 x 0r 12 09 x 0r 12 r 2 cos9 1 1 rf 2 + r? - r| / A oen 
+ + ^ &7 = r^r^l + 2 COtgl 2nr 2 ^ 

Finally 



09 1 0r 12 09 l 0r X2 09 1 0r 12 \ 90i(r x ) / r 2 cos6»! 1 r? 2 + - rf\ <9<^i(i% 
+ 7; o r 7: ^ — — ^tt, = : — 7 r - cot V\ - 



0x\ Oxi Oyi Oxi 0z\ Qx\ ) 09\ \ r i r i2 sin^x 2 2rir 2 / 09 x 

(A.36) 



The term in 



tan0i = — (A.37) 
x x 
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(v;fl( Pl ))(v 



1^12) = 




5r i2 -> 9ri2-> <9ri2 7> 

« + "7; J + — K 



dyi 



(n) 



90i 



+ 



i>i 9r 



12 



+ 



1 dr 12 \ d(j)i (r a 



9xi 9a; 1 cfyi 9xi 9zi 9x 



90! 



(A.38) 



90i 
9a; 1 
90i 
9yi 
90i 
9zi 



r\ sin 2 9 



r\ sin 2 6 1 




(A.39) 



>i9r i2 , 90i9r 12 90i 9r i2 \ 90i(r!) r 2 sin0 2 . 90i(ri) 

+ 7, q V 7, p. — ™ = '. 7T" Sin(0i - 2 ) — ^- (A.40) 



\9xi 9xi 9|/i 9xi 9zi 9a; 1 / 90i riri 2 sin0i 

The obtained expressions Eqs. (51,55) agree with the Eqs. (64-66) which had to be evaluated 
in Ref. [H] . The solution of the Eqs. (51,55) requests the use of the derivatives over spherical 
harmonics and recursion relations of cosine and sinus functions and spherical harmonics, as 
already done in Ref. 
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Figure 1: Definition and rotation of the coordinates of two electrons in an atomic center. 
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